function [allMSPs,transitionIndicies] = calculate_allMSPcomps_MASTER(data,plag2spinel,spinel2garnet,shallowP,deepP,Elements,VariablesT,VariablesX,VariablesP, garnet_coefficients, spinel_coefficients, plagioclase_coefficients)
%The process is to:
%(1) first calculate the pressure 
%(2) then using that P calculate T and X
%(3) calculate the X difference to see if the liquid is is EQ with a
%lherzolite 


%%%%%%%%%%%%%%%%%%%%

[A,TemperatureIndicies] = ismember(VariablesT, Elements);
[A,CompIndicies] = ismember(VariablesX, Elements);
TemperatureIndiciesGarnet = TemperatureIndicies; 


FittedTargetVariables = {'P' 'T' 'Qtz' 'Plag' 'Oliv' 'Cpx'}; 
ModelVariables = {'P' 'T' 'Oliv' 'Cpx' 'Plag' 'Qtz' 'P'}; 
[A,ModelVariablesIndicies] = ismember(FittedTargetVariables, ModelVariables);



data(isnan(data)) = 0; 
for i = 1:size(data,1)

%variableValues = [onesmatrix data(:,PressureIndicies)]; %add in a ones column account for the intercept



pressures = deepP:-1:shallowP;


% spinel2garnet(i)
garnetPressures =pressures(find(pressures>spinel2garnet(i)));
spinelPressures = pressures(find(pressures<=spinel2garnet(i) & pressures > plag2spinel(i)));
plagioclasePressures = pressures(find(pressures<=plag2spinel(i)));


% garnetPressures = deepP:-1:spinel2garnet(i);
% spinelPressures = spinel2garnet(i):-1:plag2spinel(i);
% plagioclasePressures = plag2spinel(i):-1:shallowP;

%Step 2 : calculate T and X

onesmatrix = ones(size(garnetPressures,2),1);
variableValues = [onesmatrix (garnetPressures)' repmat(data(i,TemperatureIndicies),size(onesmatrix,1),1)];
T = [garnetPressures' variableValues*spinel_coefficients'];
variableValues = [onesmatrix (garnetPressures)' repmat(data(i,CompIndicies),size(onesmatrix,1),1)];
X = [garnetPressures' variableValues*garnet_coefficients'];
garnetTX=[T(:,1:2), X(:,3:end)];


onesmatrix = ones(size(spinelPressures,2),1);
variableValues = [onesmatrix (spinelPressures)' repmat(data(i,TemperatureIndicies),size(onesmatrix,1),1)];
T = [spinelPressures' variableValues*spinel_coefficients'];
variableValues = [onesmatrix (spinelPressures)' repmat(data(i,CompIndicies),size(onesmatrix,1),1)];
X = [spinelPressures' variableValues*spinel_coefficients'];
spinelTX=[T(:,1:2), X(:,3:end)];


onesmatrix = ones(size(plagioclasePressures,2),1);
variableValues = [onesmatrix (plagioclasePressures)' repmat(data(i,TemperatureIndicies),size(onesmatrix,1),1)];
T = [plagioclasePressures' variableValues*plagioclase_coefficients'];
variableValues = [onesmatrix (plagioclasePressures)' repmat(data(i,CompIndicies),size(onesmatrix,1),1)];
X = [plagioclasePressures' variableValues*plagioclase_coefficients'];
plagioclaseTX=[T(:,1:2), X(:,3:end)];




singleDataPointResults = [garnetTX; spinelTX; plagioclaseTX]; 


resultsgood = singleDataPointResults(:,ModelVariablesIndicies); 
if isempty(resultsgood == 1)
    resultsgood =  NaN.*zeros(41,6);
end


allMSPs(:,:,i) =resultsgood; 



%I'm sure there must be a better way to do this, just need to keep track of
%indicies bracketing each stablity zone
transitionIndicies_temp = [1 size(garnetPressures,2) size(garnetPressures,2)+1];
transitionIndicies_temp=[transitionIndicies_temp transitionIndicies_temp(end)+size(spinelPressures,2)-1];
transitionIndicies_temp=[transitionIndicies_temp transitionIndicies_temp(end)+1  size(plagioclasePressures,2)+transitionIndicies_temp(end)];

transitionIndicies(i,:) = transitionIndicies_temp; 


% for p = 1:size(spinelPressures,2)
% variableValues = [onesmatrix spinelPressures(p)];
% spinelTX(p,:)= [spinelPressures(p) variableValues*spinel_coefficients(1:end-1,:)'];
% end
% 
% for p = 1:size(plagioclasePressures,2)
% variableValues = [onesmatrix plagioclasePressures(p)];
% plagioclaseTX(p,:)= [plagioclasePressures(p) variableValues*plagioclase_coefficients(1:end-1,:)'];
% end

%test = [garnetTX; spinelTX; plagioclaseTX];


end





%%
% modelX_normalized = bsxfun(@rdivide,modelX,nansum(modelX,2));
% dataX_normalized = bsxfun(@rdivide,dataX,nansum(dataX,2));
% NRMSD= sqrt(sum((modelX_normalized-dataX_normalized).^2, 2));
%  fits = [RMSD,NRMSD] ;
% figure(); hold on; plot(NRMSD,RMSD,'*k')
% xlabel('NRSMD')
% ylabel('RMSD')

%%
%mat2clip(dataX)



end

